classdef Element < hgsetget
    %Clase que representa un elemento de la malla
    %Hereda del handle hgsetget para poder utilizar los métodos set y get
        
    properties
        icone
        K_diffusion %matriz de rigidez (incluye cb robin si la hay)
        b_diffusion %flujo difusivo (incluye cb robin si la hay)
        K_convection %matriz convectiva
        b_source %fuente
        K_reaction %matriz reactiva
        K_temporal %matriz de masa temporal
        b_temporal %carga temporal
        q_faces %nodos(1d) o caras(2d) del borde con condicion de flujo
        h_faces %nodos(1d) o caras(2d) del borde con condicion mixta
    end
    
    methods
        function [obj] = Element(obj)
            
        end
        
        function [obj] = assemble(obj)
        end
        
        function [obj] = assemble_diffusion(obj)
            global data;
            % de cada nodo que conforma el elemento, tomo las coordenadas X
            % e Y y su difusividad
            xx = cell2mat(get(data.mesh.nodes(obj.icone),{'x'}));
            k_ele = cell2mat(get(data.mesh.nodes(obj.icone),{'k'}));
            
            %numero de nodos en el elemento y numero de dimensiones del
            %problema
            [nen,ndm] = size(xx);
            
            [pg,wpg] = data.sF.getGaussPointsAndWeights();
            npg = size(pg,1);
            
            % Ensamblo la matriz elemental
            obj.K_diffusion = zeros(nen,nen);
            obj.b_diffusion = zeros(nen,1);
            
            %loop sobre cada punto de gauss
            for ipg=1:npg
                %mapeo las coordenadas al elemento master correspondiente y traigo
                %los gradientes       
                [dN,N,J0,detJ] = data.sF.getShapeFunctions(data.param.Gtype_shfunc,xx,pg(ipg,:));
                
                %propongo "k" constante por elemento:
                k = mean(k_ele);
                
                for inod=1:nen
                    for jnod=1:nen
                        for idm=1:ndm
                            for jdm=1:ndm
                                % Aij = -k*int(gradNi*gradNj)
                                obj.K_diffusion(inod,jnod) = obj.K_diffusion(inod,jnod) - k * (dN(inod,jdm)*J0(jdm,idm) * dN(jnod,jdm)*J0(jdm,idm)) * wpg(ipg)/detJ;
                            end
                        end
%                         if ndm==2  
%                             A(i,j) = -k'*[ (dN(i,:)*J0(:,1))*(dN(j,:)*J0(:,1))'; (dN(i,:)*J0(:,2))*(dN(j,:)*J0(:,2))'] * wpg(ipg)/abs(detJ);
%                         else
%                             A(i,j) = -k * (dN(i)*J0*dN(j)*J0) * wpg(ipg)/abs(detJ);
%                         end
                    end
                end
            end

            nfaces = nen; %para todo elemento, tengo el mismo numero de nodos que de caras...?
            if ndm==2
                if nen==3 %elemento triangular lineal -> lineal en el borde
                    faces = [1 2; 2 3; 3 1];
                    nen_face=2;
                elseif nen==6  %elemento triangular cuadrático -> cuadratico en el borde
                    faces = [1 2 3; 3 4 5; 5 6 1];
                    nen_face=3;
                elseif nen==4;  %elemento rectangular lineal  -> lineal en el borde
                    faces = [1 2; 2 3; 3 4; 4 1];
                    nen_face=2;
                end
            elseif ndm==1
                if nen==2 %elemento lineal ->los dos nodos pertenecen al borde del elemento
                    faces = [1;2];
                elseif nen==3  %elemento cuadratico -> el primero y el ultimo nodo pertenecen al borde del elemento
                    faces = [1;3];
                end
                nen_face = 1;
            end
            
            for iface=1:nfaces %para cada nodo o cara reviso si debo aplicar CB
                
                %primero verifico condicion Neumann
                if find(iface==obj.q_faces) %el nodo o la cara tiene condicion Neumann de flujo
                    q=0;
                    for I=1:nen_face
                        q = q + get(data.mesh.nodes(obj.icone(faces(iface,I)),'q'));
                    end
                    q = q/nen_face; %valor de q promediado en la cara
                
                    %la longitud de la cara es la distancia entre el
                    %primer y ultimo nodo
                    if ndm==2
                        l = sqrt((xx(faces(iface,1),1)-xx(faces(iface,end),1))^2 + (xx(faces(iface,1),2)-xx(faces(iface,end),2))^2);
                    elseif ndm==1
                        l = 1;
                    end
                    
                    for I=1:nen_face
                        obj.b_diffusion(I) = obj.b_diffusion(I) + q*l/nen_face;
                    end
                end

                if find(iface==obj.h_faces) %el nodo o la cara tiene condicion Neumann de flujo
                    h=0; phi_inf=0;
                    %utilizando
                    for I=1:nen_face
                        h = h + get(data.mesh.nodes(obj.icone(faces(iface,I)),'h'));
                        phi_inf = phi_inf + get(data.mesh.nodes(obj.icone(faces(iface,I)),'phi_inf'));
                    end
                    h = h/nen_face; %valor de h promediado en la cara
                    phi_inf = phi_inf/nen_face;  %valor de phi_inf promediado en la cara
                    
                    if ndm==2
                        %la longitud de la cara es la distancia entre el
                        %primer y ultimo nodo
                        l = sqrt((xx(faces(iface,1),1)-xx(faces(iface,end),1))^2 + (xx(faces(iface,1),2)-xx(faces(iface,end),2))^2);
                        %aporte del termino mixto a la matriz de difusion,
                        %integrar sobre el contorno
                        for I=1:nen_face
                            obj.b_diffusion(faces(iface,I)) = obj.b_diffusion(faces(iface,I)) - h*phi_inf*l/nen_face;
                            for J=1:nen_face
                                inod = faces(iface,I);
                                jnod = faces(iface,J);
                                %integro con puntos de gauss en una dimension
                                %menor (2d->1d)
                                %xy->t parametrizo las dos dimensiones con
                                %un parametro t que vaya desde 0 a l
                                %(t=[0,l]) e integro allí
                                t = [0:l/(nen_face-1):l]';
                                [pg,wpg] = data.sF.getGaussPointsAndWeights(ndm-1);
                                for ipg=1:npg
                                    [dN,N,J0,detJ] = data.sF.getShapeFunctions(data.param.type_Gshfunc,t,pg(ipg,:));
                                    obj.K_diffusion(inod,jnod) = obj.K_diffusion(inod,jnod) - h * N(inod) * N(jnod) * wpg(ipg)/abs(detJ);
                                end
                            end
                        end
                    elseif ndm==1
                        l = 1;
                        obj.K_diffusion(inod,jnod) = obj.K_diffusion(inod,jnod) - h*l/nen_face;
                    end
                    
                end

            end
        end
        
        function [obj] = assemble_convection(obj)
            %analizar que funciones de forma utilizar... galerkin standart
            %streamline diffusion, etc
        end
        
        function [obj] = assemble_reaction(obj)
            global data;
            % de cada nodo que conforma el elemento, tomo las coordenadas X
            % e Y y su constante de reaccion
            xx = cell2mat(get(data.mesh.nodes(obj.icone),{'x'}));
            r_ele = cell2mat(get(data.mesh.nodes(obj.icone),{'r'}));
            
            %numero de nodos en el elemento y numero de dimensiones del
            %problema
            [nen,ndm] = size(xx);
            
            [pg,wpg] = data.sF.getGaussPointsAndWeights();
            npg = size(pg,1);
            % Ensamblo la matriz elemental
            obj.K_reaction = zeros(nen,nen);
            
            %loop sobre cada punto de gauss
            for ipg=1:npg
                %mapeo las coordenadas al elemento master correspondiente y traigo
                %los gradientes       
                [dN,N,J0,detJ] = data.sF.getShapeFunctions(data.param.Gtype_shfunc,xx,pg(ipg,:));
                
                %propongo "r" constante por elemento:
                r = mean(r_ele);
                
                for inod=1:nen
                    for jnod=1:nen
                        for idm=1:ndm
                            for jdm=1:ndm
                                % Aij = int(r*Ni*Nj)
                                obj.K_reaction(inod,jnod) = obj.K_reaction(inod,jnod) - r_ele(inod) * N(inod) * N(jnod) * wpg(ipg)/abs(detJ);
                            end
                        end
                    end
                end
            end
        end
        
        function [obj] = assemble_source(obj)
            global data;
            % de cada nodo que conforma el elemento, tomo las coordenadas X
            % e Y y su fuente
            xx = cell2mat(get(data.mesh.nodes(obj.icone),{'x'}));
            s_ele = cell2mat(get(data.mesh.nodes(obj.icone),{'s'}));
            
            %numero de nodos en el elemento y numero de dimensiones del
            %problema
            [nen,ndm] = size(xx);
            
            [pg,wpg] = data.sF.getGaussPointsAndWeights();
            npg = size(pg,1);
            
            % Ensamblo el residuo elemental
            obj.b_source = zeros(nen,1);
                
            %loop sobre cada punto de gauss
            for ipg=1:npg
                %mapeo las coordenadas al elemento master correspondiente y traigo
                %los gradientes       
                [dN,N,J0,detJ] = data.sF.getShapeFunctions(data.param.Gtype_shfunc,xx,pg(ipg,:));
                
                
                %propongo "s (fuente)" constante por elemento:
                %s = mean(s_ele);

                for inod=1:nen
                    obj.b_source(inod) = obj.b_source(inod) + s_ele(inod) * N(inod) * wpg(ipg)/detJ;
                end
            end
        end
        
        function [obj] = assemble_temporal(obj)
        end
        
    end
    
end

